Dynamical law of the phase interface motion in the presence of crystals nucleation

In this paper, we develop a theory of solid/liquid phase interface motion into an undercooled melt in the presence of nucleation and growth of crystals. A set of integrodifferential kinetic, heat and mass transfer equations is analytically solved in the two-phase and liquid layers divided by the moving phase transition interface. To do this, we have used the saddle-point method to evaluate a Laplace-type integral and the small parameter method to find the law of phase interface motion. The main result is that the phase interface Z propagates into an undercooled melt with time t as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z(t)=\sigma \sqrt{t}+\varepsilon \chi t^{7/2}$$\end{document}Z(t)=σt+εχt7/2 with allowance for crystal nucleation. The effect of nucleation is in the second contribution, which is proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{7/2}$$\end{document}t7/2 whereas the first term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim \sqrt{t}$$\end{document}∼t represents the well-known self-similar solution. The nucleation and crystal growth processes are responsible for the emission of latent crystallization heat, which reduces the melt undercooling and constricts the two-phase layer thickness (parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi <0$$\end{document}χ<0).

Growth rate constant, -

Abbreviations
Heav Heaviside function WVFZ Weber-Volmer-Frenkel-Zel' dovich nucleation kinetics The dynamics of crystallization fronts propagation have attracted the attention of researchers for more than 130 years, starting with Stefan's famous works on the freezing of water with a flat front 1,2 . Nowadays, problems with a moving boundary separating different phases of the aggregate state of matter bear his name. The rich variety of nonlinear dynamics of interfacial boundaries has attracted essential interest in applied mathematics, which deals with various moving boundary problems in geophysics, materials science, nonlinear physics, and heat and mass transfer theory (see, among others, [3][4][5][6].
It is well-known that the phase transition boundary in Stefan's problem moves according to the self-similar law (it is proportional to the square root of time). Such a law is preserved even in the presence of an extended phase transition region separating purely solid and liquid phases in binary 7 and three-component 8,9 melts. For example, in the case of binary melts, there is one two-phase region and two interphase boundaries moving in proportion to the square root of time t 7 . In the case of ternary melts, there are two two-phase regions (main and cotectic) and three interphase boundaries, moving according to the law ∼ √ t with various parabolic growth rate constants 9 . Note that self-similar solutions take place only in the semi-infinite region whose boundary temperature is kept constant. At a given heat flow at this boundary or its active cooling by air or water flow, there is no self-similarity of the phase interface motion 10 . There is also no self-similarity of motion of the interfacial boundary (boundaries) when considering nonlinear heat conduction and impurity diffusion equations, heat and mass sources in these equations, heat exchange with the environment, melt convection, kinetics, stochastic fluctuations as well as other factors (see, among others [11][12][13][14][15][16]. In this study, we derive the law for phase interface motion with allowance for nucleation and growth of crystals in front of it. These effects significantly change the previously known self-similar law due to the thinning of the two-phase layer as a result of additional latent heat released by the growing crystals. The theory developed can be used in describing the propogation of various solid/liquid interfaces in metastable media and experimental data on Al-rich Al-Ni alloys showing an anomalous solidification behaviour: the solid-liquid interface velocity slows down as the undercooling increases 17,18 .

The model
First of all, let's assume a one-component melt filling a semi-space η > 0 ( Fig. 1). We will also assume that the initial temperature T l of the melt exceeds its crystallization temperature T p at the time τ = 0 . The temperature T = T 0 at the solid wall η = 0 then jumps to a value of T < T p . This generates the melt undercooling T = T p − T propagating into the melt phase 0 < η < �(τ) . Here, �(τ ) stands for the interface boundary with temperature T = T l = T p (T and T l denote the temperatures in the moving domains 0 < η < �(τ) and η > �(τ ) ). The undercooled domain 0 < η < �(τ) is filled with nucleating particles whose nucleation rate is I(�T) ( T = 0 at the interface η = �(τ ) ). As a result of crystal growth, the melt undercooling is partially reduced due to the release of latent heat of crystallisation. This phenomenon is the reason for slowing down the interfacial boundary movement into the liquid melt.
For the sake of simplicity, we consider the growth of an ensemble of spherical particles whose distribution function φ(η, r, τ ) satisfies the first-order kinetic equation. In addition, we have the thermal conductivity equations in both the domains with allowance for the heat source in the two-phase region as a result of crystal growth. The governing equations read as www.nature.com/scientificreports/ Here Q v is the latent heat of phase transition, ρ and ρ l are the densities of two-phase and melt phases, c and c l are their heat capacities, and and l are their heat conductivity coefficients, respectively (subscript l denotes liquid region).
These equations should be supplemented with the conditions at r = 0 , τ = 0 , η = 0 , η → ∞ , and the phasetransition boundary η = �(τ ) , of the form of Let us especially underline that the first expression (4) gives the flux of nucleating crystals that appear within the two-phase region.
where β k and n are constants. These parameters can be found from experimental data on crystal growth in a certain undercooled liquid. For the sake of simplicity, we consider the case when the nucleation rate depends only on the melt undercooling. Dealing the case of the Weber-Volmer-Frenkel-Zel' dovich (WVFZ) nucleation kinetics, we have 23 where T 0 represents the initial melt undercooling. Another frequently used kinetics is Meirs power law [24][25][26] It is significant that the constant parameters I * and p entering expressions (8) and (9) are different.

Analytical solutions to moving-boundary problem
The aforementioned model represents a set of integrodifferential equations, initial and boundary conditions with a moving phase interface �(τ ) . Below we develop a method to its analytical solution, which is based on the saddle-point technique to evaluate a Laplace-type integral 27,28 and the small parameter technique to solve a nonlinear heat transfer equation for the melt undercooling.
To simplify the matter, we use dimensionless parameters as follows www.nature.com/scientificreports/ where T 0 = T p − T 0 and T l = T p − T l represent the initial and current melt undercoolings. The dimensionless distribution function F of growing crystals in a two-phase region is defined by the following problem (this problem can be easily obtained from Eqs. (1), (4), (7)-(9) rewritten in dimensionless form) where Applying the Laplace integral transform to the problem (11) and (12) with respect to t, we find the particleradius distribution function of the form where and Heav(·) is the Heaviside function. Now substituting dimensionless variables and parameters (10) into Eqs. (2), (3) and boundary conditions (5) and (6), we obtain where w > 0 in the two-phase layer, and w l < 0 in the liquid layer. Also, for the sake of simplicity, we assume that = l , ρ = ρ l , and c = c l . Let us especially emphasize that the new variable ν has been used in Eq. (15) so that x(z, ν) = x(z, t) − s , and To evaluate the integral in the r.h.s. of Eq. (15) we apply the saddle-point technique for a Laplace-type integral 27,28 . The ν-derivatives of the function g are negative for both nucleation kinetics under consideration since the melt undercooling w is a decreasing function of time: ∂g/∂ν = (dg/dw)∂w/∂ν < 0 . It means that the maximum point of the function g lies at the left boundary at ν = 0 . To calculate the first non-zero derivative of g(z, ν) with respect to ν we use the same Eq. (15). It is not difficult to show that the first three derivatives are zero and only the fourth one is non-zero for both nucleation kinetics. Namely, g (4) = −12b (WVFZ) and g (4) = −6b (Meirs) at ν = 0 . Now keeping in mind only the main contribution in the asymptotic expansion of the integral (15), we come to Ref. 27,28 (10) WVFZ nucleation mechanism ln w, Meirs nucleation mechanism. www.nature.com/scientificreports/ where ε = A −n ; µ = 2 and µ = 1 for the WVFZ and Meirs mechanisms. Our estimations show that ε ≪ 1 for typical metallic melts 29,30 . This fact enables us seeking for the solution of Eq. (19) expanding the rescaled undercooling u in a series in small parameter ε as Substituting (22) into (19) and (20), expanding conditions at z = Z(t) in series and equating terms with the same power of ε , we arrive at the following form of solutions where σ and χ represent the constant parameters characterizing the interface position Z(t). The functions � 0 (ζ ) and � 1 (ζ ) satisfy the following equations and boundary conditions where The analytical solution to the model (24)-(26) takes the form Here, parameters σ and χ satisfy the equations Let us especially note that expressions (27) and (28) define the melt undercooling w > 0 within the two-phase region at 0 < z < Z(t) and the undercooling w l < 0 in pure melt at z > Z(t).

Behaviour of solutions
Our analytical solutions are illustrated in Figs. 2, 3 and 4 for parameters typical for metallic melts 30 . First of all, phase interface dynamics, shown in Fig. 2 with allowance for nucleation and growth of crystals, essentially differ from the case without particles in a two-phase layer. This purely frontal case demonstrated by the dotted curve in Fig. 2 is described by the law of the square root of time, i.e. Z 0 (t) ∼ √ t. Such dynamics are the property of so-called self-similar crystallization processes (see, among others 7,31 ). Nucleation and growth of solid particles within the two-phase layer change this dynamical law drastically. From (20) ∂u l ∂t =γ ∂ 2 u l ∂z 2 , z > Z(t), t > 0, www.nature.com/scientificreports/ a certain point in time, the law of motion of the phase interface becomes a decreasing function of time. This is caused by the fact that growing crystals produce the latent heat of phase transition, which partially reduces the undercooling in a two-phase layer and, thus, its thickness. In addition, the power of nucleation rate n significantly affects the movement of the phase transition interface Z(t) (compare the solid and dashed curves in Fig. 2). The greater n (higher crystal growth rate according to expression (7)) the greater the two-phase layer thickness. Moreover, increase of n shifts the maximum point of phase transition boundary Z(t) towards higher values of time t, i.e. the compression of the two-phase layer occurs later with increasing n. Figure 3 illustrates the melt undercooling in the two-phase ( w > 0 ) and liquid ( w l < 0 ) layers. As is easily seen, nucleation and crystal growth processes substantially constrict a two-phase layer and accelerate its desupercooling dynamics when time t increases (compare the green dotted and blue solid curves in Fig. 3 shown for different time instants t). In addition, the melt undercooling in liquid becomes lower with increasing ζ and t. As this takes place, the first correction w 1 = u 1 /A to the main contribution w 0 = u 0 /A substantially influences  www.nature.com/scientificreports/ the desupercooling dynamics as compared with w 0 . Such a dynamics could be the key to explain the anomalous U-shape behaviour of the "recalescence front velocity-melt undercooling" curve in Al-rich Al-Ni alloys 17,18 . The particle-radius distribution function (13) at different points z in the two-phase layer is shown in Fig. 4 at a certain point in time. As is easily seen, this function represents a bell-shaped curve decreasing its amplitude with increasing the spatial coordinate z in a two-phase region (when approaching the two-phase layer-liquid phase boundary). This is due to the fact that the melt undercooling increases with decreasing z (when approaching the two-phase layer boundary z = 0 ). Such a bell-shaped behaviour is in agreement with typical particleradius distribution in undercooled and supersaturated liquids (see, among others, recently published review on nucleation 32 ).

Conclusion
In summary, we develop a theory of solid/liquid phase interface motion in the presence of nucleation and particle growth processes in an undercooled layer. This layer together with evolving crystals, which partially compensate for the undercooling, propagates into pure melt with the velocity depending on crystal growth and nucleation rates. The main result of our study is the dynamical law of the phase interface boundary motion Z(t) = σ √ t + εχt 7/2 (two coefficients σ and χ are found analytically). This law substantially differs from the case without nucleation and growth of crystals, which defines the phase interface as Z 0 (t) = σ √ t . The negative coefficient χ leads to the narrowing of an undercooled layer starting from a certain point in time. This is caused by the effect of latent heat emission, which reduces the melt undercooling and constricts its spatial size (the two-phase layer thickness).
The present theory can be extended to take into account fluctuations in the growth rates of individual particles in an undercooled two-phase layer. This effect raises the order of the kinetic equation in the spatial variable (see, among others 32,33 ). Such consideration can modify the law of motion of the phase interface. However, the slowing down of its motion detected in this paper in comparison with the self-similar case ( Z 0 (t) = σ √ t ) should be preserved, since it is caused by the compensation of undercooling due to the release of the latent crystallization heat.

Data availability
All data generated or analysed during this study are included in this published article.